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A correlated quantum many-body method is applied to describe resonance states of atomic Bose- 
Einstein condensates (BEC) in a realistic shallow trap (as opposed to infinite traps commonly used). 
The realistic van der Waals interaction is adopted as the interatomic interaction. We calculate 
experimentally measurable decay rates of the lowest quasi-bound state in the shallow trap. The most 
striking result is the observation of a new metastable branch besides the usual one for attractive 
BEC in a pure harmonic trap. As the particle number increases the new metastable branch appears, 
then gradually disappears and finally usual metastable branch (associated with the attractive BEC 
in a harmonic trap) appears, eventually leading to the collapse of the condensate. 
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I. INTRODUCTION 



In experiments with Bose-Einstein condensation, the 
evolution of atomic cloud and its instability strongly 
depends on the external confinement, which is usually 
chosen as either isotropic or anisotropic pure (i.e. of 
infinite extent) harmonic potential. But, in the actual 
experimental setup, the trap is of finite extent. During 
the last few years, attention has been shifted to shallow 
optical dipole traps [l|. As the quadratic trapping 
potential takes the shallow Gaussian envelope form, 
the anharmonicity of the potential must be taken into 
account. BEC in such a shallow trap of finite width 
supports resonance states which are quasi-bound. In 
such an experimental trap, the decay mechanism of 
the condensate becomes an important issue, as the 
condensate can escape from the trapping potential by 
quantum tunneling through intermediate barriers, in 
addition to the usual collapse of attractive condensates. 
Several attempts have been made to calculate life 
time of quasi-bound states and to study the transition 
from a resonance to a bound state by solving the 
Gross-Pitaevskii equation (GPE) using contact delta 
interaction The GPE is based on the mean- field 

description and ignores correlations in the many-body 
wave function. But near the criticality, the energy of the 
resonance state becomes very close to the barrier height 
and the condensate becomes highly correlated. Thus 
incorporation of interatomic correlation becomes impor- 
tant. Naturally the full quantum many-body treatment 
incorporating a realistic interatomic interaction is indeed 
necessary. In particular the present experiments consider 
only a finite number of atoms in the external trap. Such 
condensates are quantum depleted, which again deserve 
a full quantum many-body treatment. The external 
shallow potential can be modeled by a quadratic plus a 
quartic potential, viz., V(r) = ^mui 2 r 2 + Ar 4 , where A 
is the anharmonic parameter. In our earlier work in this 
direction, we investigated the ground state properties 



in such a trap. We also observed dramatic change in 
the stability factor Nc ^ asc ^ (where N cr is the critical 
number of atoms beyoncf which the attractive condensate 
collapses) due to the anharmonicity Q- 



Thus the aim of our present work is to employ a 
correlated quantum many-body approach incorporating 
a finite range realistic interatomic interaction to study 
resonance states in shallow traps. We calculate decay 
rates of quasi-bound condensates and study how the 
interatomic interaction and the anharmonicity interfere 
with the decay process. Deviations from earlier works 
which used the mean-field equation can be attributed 
to two-body correlations and finite range of the realistic 
interaction. Due to the realistic nature of our calculation 
we expect results which are relevant to experiments. 
We observe resonance states associated with an unusual 
metastability for the repulsive BEC. The most striking 
result is the observation of two metastable branches 
for an attractive BEC. The attractive BEC containing 
N atoms in a pure harmonic trap is associated with 
a metastable branch, which ultimately collapses for 
N > N cr . In the present study we observe that with 
increase in particle number, a new metastability first 
appears, then gradually disappears and finally the usual 
metastability (associated with attractive BEC in a 
pure harmonic trap) appears, eventually leading to the 
collapse. We also study macroscopic quantum tunneling 
and calculate decay rates for these two branches. 



In Sec. II we introduce the methodology which contains 
the many-body approach used in this work based on the 
correlated potential harmonic expansion method. Sec. Ill 
contains numerical results and discussions. Conclusions 
are drawn following a summary of our work in Sec. IV. 
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II. CORRELATED POTENTIAL HARMONIC 
EXPANSION METHOD (CPHEM) 

We adopt the potential harmonic expansion method 
with a short range correlation function (CPHEM) which 
has already been established as a very successful and 
useful technique for the study of dilute BEC 0413] ■ 
Here we describe the technique briefly. 

We consider that a system of A = (N + 1) identical 
bosons, each of mass m, is confined in an external trap 
(V trap (r)) which is modeled as a harmonic potential (of 
frequency uj) perturbed by a quartic term. The time 
independent quantum many-body Schrodinger equation 
is written as: 



2m 



A 



? Vtrap 



i=l »=1 

-E]^(x u ...,x A ) =0 



i,j>i 



(1) 



Where E is the total energy of the system, V(xi — Xj) is 
the two-body potential and Xi is the position vector of 
the ith particle. We define a set of N Jacobi vectors as: 



and are inert spectators. Therefore, only two-body cor- 
relations in the many-body wave function and pair-wise 
interactions among atoms are important. For the (ij)- 
interacting pair, we define a hyperradius OjA for the re- 
maining (N — 1) noninteracting bosons as[12j 



Pij 



\ 
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E 
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C 2 



(5) 



so that r 2 = 



A hyperangle <p is introduced 



through Xij = rcoscp and = rsin0. The full set of 
3iV hyperspherical variables are chosen as: 

a) r, </>, ip ($ and ip are the polar angles of fy corre- 
sponding to the interacting pair) and, 

b) (37V — 4) hypcrangular variables associated with the 
remaining (N — 1) inert spectators. Out of these, 2(N — 
1) are the polar angles associated with (N — 1) Jacobi 
vectors and another (N — 2) 'hyperangles' define their 
relative lengths. 

In the Potential Harmonic Expansion 
Method (PHEM), we decompose the many-body 
wave function^) into Faddeev components (pij for the 
(«j)-interacting pair as 



— — -{Xi+1 - - > x 
i + 1 i 



) (i = l,...N)- (2) 



The center of mass coordinate is R 



i=l 



As the 



labelling of particles is arbitrary, we choose the relative 
separation of (r/)-interacting pair {xij = Xi — Xj) as Cat. 
We define the hyperradius r as 



N 



Ec 2 



a E x l 

i,j>i 
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(3) 



where fi = Xi — R is the position vector of z-th particle 
with respect to the center of mass of the system. The 
relative motion of the bosons is given by 



/v 



5~*V?. +Vtrav + Vint{(l, Ov) 



-E R *(Ci,...,Cat) = 



(4) 



Vtrap is the effective external trapping potential and Vi n t 
is the sum of all pair- wise interactions. Er is the rel- 
ative energy of the system i.e. E = En + ^hu>. The 
laboratory BEC is designed to be very dilute, so that 
the probability of three or more atoms to come within 
the range of interaction is negligible. This is done in 
laboratory experiments, in order to avoid molecule for- 
mation through three body collisions. Hence, when the 
(ij) pair interacts, the rest of the atoms are far apart 



* = E 



(0) 



i,j>i 



Then Eq. (4) can be expressed as 



[T + Vtrap - E R ] (pij = ~V(Xij) 2J <t>kl 



(7) 



kl>k 



N 



where T = — ^"E/^Ci' ^ mcc only two-body correlations 

i=l 

are relevant, the Faddeev component <f>ij is independent 
of the coordinates of all the particles other than the in- 
teracting pair. Hence the total hypcrangular momentum 
quantum number as also the orbital angular momentum 
of the whole system (comprising of all bosons) is con- 
tributed by the interacting pair only. We expand (pij in 
the subset of hyperspherical harmonics (HH) necessary 
for the expansion of V{xij). Thus 



4>ij { x ij i r , 



^E^+^^w- (8) 



K 



where denotes the full set of hyperangles in the 3./V- 
dimensional space corresponding to the (ij)-interacting 
pair. This new basis set {T^x+z^tv)} is called po- 
tential harmonics (PH) basis and it is independent of 
(Ci, C/v-i)- Thus the total angular momentum of the 
condensate and its projection are simply I and m, which 
comes from the interacting pair and all other quantum 
numbers coming from (-/V— 1) non- interacting bosons are 
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kept frozen. An analytic form of the potential harmon- 
ics can be found in Ref [12]. Substitution of Eq. (8) into 
Eq. (7) and taking projection on a particular PH gives Q 



" + + Vtra P (r) E R ] U< K (r) 

^2fhiVKK>(r)U l K ,(r) = 0- 



K> 



when Vkk' is the potential matrix and is given by 



(9) 



V KK >(r) 



where Ck — 2K + I 



(10) 

3N 2 3 and ff Cl is the overlap of 
the PH for ij partition with the sum of PHs of all par- 
titions Q. K is the grand orbital quantum number in 
3iV dimensional space. All other intermediate angular 
momentum quantum numbers take zero eigen values. As 
the number of active variables is now only four (global 
hyperradius r and three other for Xij ) for any number of 
bosons, the numerical complication is greatly simplified. 
Eq. (9) can be put in symmetric form as 



^ + V trap (r) + ^{C(£ + 1) 



+4K(K + a + /3 + l)}-E R 
^2fKiV K K'(r)fK'iU K 'i(r) = 



U m {r) 
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where C = I + ^=fi, U Kl 
(3 = 1 + 1/2. 



fmu l K {r). a 
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and 



In experimentally achieved BEC, as the energy of the 
interacting pair is extremely small, the two-body interac- 
tion is reproduced by the s-wave scattering length (a sc ). 
Positive (negative) a sc corresponds to a repulsive (at- 
tractive) condensate. In the Gross-Pitaevskii equation, 
the interatomic interaction is a contact interaction and 
is absolutely determined by its strength a sc only. Thus 
the two-body potential is purely repulsive or purely at- 
tractive, depending on the sign of a sc . But van der Waals 
potential has two terms, one part represents a strong re- 
pulsion (usually represented by a hard core of radius r c ) 
at very short separation and the other part goes as — -^f 

and a sc can be either positive or negative depending on 
the value of r c Q. Thus we determine a sc by solving 
the zero-energy two-body Schrodinger equation for the 
two-body wave function rj(xij) 



h 2 l 



drj(xj 



+ V(x ij )r ) (x ij )=0- (12) 



The value of a sc is determined from the asymptotic part 
of rj(xij) [li|. The zero-energy two-body wave function 
r)(xij) is also a good representation of the short range 
behavior of <j>ij and is taken as the two-body correlation 
function in the PH expansion basis, to improve the rate 



of convergence [l4[. Thus in the CPHEM, we replace 
Eq. (8) by 



>ij{xij,r) = r ^ , 



£P&*(n£)«Ur)»7(stf) 



K 



(13) 

Introduction of rj{xij) enhances the rate of convergence 
of the PH expansion dramatically. This has been ac- 
tually verified in our numerical calculation. In our nu- 
merical procedure we solve Eq. (12) for the zero energy 
two-body wavefunction rj(xij) in the chosen two-body po- 
tential V(xij). We adjust the hard core radius r c , such 
that a sc has the desired value 0, [l3[. This r](xij) is then 
used in Eq. (13) and the potential matrix becomes 




V KK >(r) = (h K 

S+l{pf{z)v( 



Wi(z)}dz ■ 
(14) 

where (z) is the Jacobi polynomial, and its norm and 
weight function are h^f and Wi(z) respectively [15j . We 
truncate the if-sum in Rq. (13) to an upper limit K max 
providing the desired convergence. Finally the coupled 
differential equation (CDE), Eq. (11), is solved by the 
hyperspherical adiabatic approximation (HAA) [16|. In 
HAA, one assumes that the hyperradial motion is slow 
compared to the hyperangular motion. Hence the lat- 
ter is separated adiabatically and solved for a particular 
value of r, by diagonalizing the potential matrix together 
with the diagonal hypercentrifugal repulsion in Eq. (11). 
The lowest eigenvalue, uio(r) (corresponding eigen col- 
umn vector being XKo( r ))i provides the effective poten- 
tial for the hyperradial motion. We choose the lowest 
eigen potential (uJo(r)) as the effective potential in which 
the entire condensate moves as a single entity. Thus in 
HAA, the approximate solution (the energy and wave 
function) of the condensate is obtained by solving a sin- 
gle uncoupled differential equation 



h 2 d 2 

+ - E R 



Co(r) = 



(15) 



subject to appropriate boundary conditions on C,q{t). 
The function CoM is the collective wave function of the 
condensate in the hyperradial space. The lowest lying 
state in the effective potential u!q(t) corresponds to the 
ground state of the condensate. The total energy of 
the condensate is obtained by adding the energy of the 
center of mass motion (^hu>) to E R , 

The main advantages of our CPHEM are: 

i) Potential harmonic basis keeps all possible two-body 
correlations and yet the number of variables is reduced 
to only four for any number of bosons in the trap. So 
in spite of incorporating all the two-body correlations, 
we can treat quite a large number of atoms in the trap 
without much numerical complication. 

ii) We can use a realistic interatomic interaction like 
the van der Waals potential having a finite range, which 
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itself takes care of the short range repulsion and inter- 
atomic correlations. 

iii) Unlike the GP equation, CPHEM does not have any 
pathological singularity, since the two-body interaction is 
a realistic one and has a strong short-range repulsion. 

Thus the CPHEM reveals the realistic picture. Clearly 
it is an improvement over the GP equation. Finally, 
by using the HAA, we reduce the multi-dimensional 
problem into an effective one-dimensional one in hy- 
perradial space and the effective potential uo(r) of this 
one-dimensional problem provides a clear qualitative 
picture and a quantitative description of the system. 



quartic confinement becomes more tight, the frequency 
will increase for repulsive BEC and the attractive 
BEC will again be associated with a metastability Q. 
These have been studied earlier both in one and three 
dimensions [22|, [23j] . However the most dramatic features 
are expected for A<0, and the potential can be easily 
approximated as V (r) = ^r 2 exp(— cr 2 ) with A ~ |. 
Our present calculation will consider only A<0. 



A. Repulsive BEC 



III. RESULTS 

We choose the interatomic potential as the van der 
Waals potential with a hard core of radius r c , viz. 
V(xij)= oo for Xij < r c and = ——§- for Xij > r c . 

The strong short range repulsion is parameterized by 
the hard core and the strength (Cg) is known for a 
given type of atom, e.g., C§ = 6.4898 x 10~ n o.u. for 
Rb atoms [l3|. The value of r c is adjusted to get the 
desired value of a sc . In oscillator unit (o.u.) length 

and energy are given in the units of a^ Q = and 
htu respectively. As C 6 -> 0, the potential becomes a 
hard core potential and r c coincides with the s-wave 
scattering length. As detailed in the previous section, 
we solve the zero-energy two-body Schrodinger equation 
for the interacting pair to get the value of r c , which 
corresponds to the experimental scattering length a sc . 
With a tiny change in r Cj a sc may change by a large 
amount, including sign [13j . With each additional 
change in sign, the potential supports an extra two-body 
bound state resulting in an additional node in rj{xij). 
Thus the choice of r c is very crucial. We choose r c such 
that it corresponds to the zero-node in the two-body 
wave function for attractive interaction and one node 
for repulsive interaction (lp| . For a repulsive BEC we 
choose 87 Rb atoms with a sc = .00433 o.u. as in the 
JILA trap [l?}- For an attractive BEC, we choose 
a sc = —1.832 x 10~ 4 o.u., which is one of the values as 
reported in the controlled collapse experiment of 85 Rb 
atoms [ill, [III. In both cases r c is determined by the 
method discussed above. 

In the optical dipole trap, the trapping potential takes 
the shallow Gaussian form and the external trap is given 



by V(r) 



1^,2 



+ Xr . For A > the frequency is blue 



shifted and for A < the frequency is red shifted. In 
the experiment [20L l2lj . quartic confinement is created 
with a blue detuned Gaussian laser directed along the 
axial direction. The nonrotating condensate was cigar 
shaped and the strength of the quartic confinement 
was w 10 -3 . In the present study we choose A as a 
controllable parameter and |A| <C 1. For A>0, as the 



For harmonic trapping with repulsive interaction, the 
condensate is always stable for any number of bosons. 
However due to the presence of anharmonicity we ob- 
serve a new and different metastablity of the condensate. 
In Fig. 1 we plot the effective potential wo(r) as a func- 
tion of r for 500 atoms of 87 Rb in a shallow trap corre- 
sponding to A = —2 x 10~ 5 o.u. and a sc = .00433 o.u.. 
We observe a dramatic change in the effective potential 
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FIG. 1: (color online) Plot of the effective potential coo(r) 
against r (both expressed in appropriate o.u.) for N = 500 
atoms of 87 Rb with a sc = .00433 o.u. and A = -2 x 10~ 5 o.u. 
The metastable region (MSR) and the intermediate barrier 
on the right side (RIB) are indicated in the figure. 



from that of a purely harmonic trap: a metastable re- 
gion (MSR) appears bounded by an intermediate barrier 
on the right side (RIB), beyond which wo(r) decreases 
gradually, where the quartic term dominates over the 
quadratic term. In our earlier work Q, we studied how 
the ground state properties and the low energy collective 
excitations get modified due to anharmonicity and calcu- 
lated the stability factor JVe ^ asc in such a shallow trap. 
In the present work, we calculate the decay rate of quasi- 
bound states in the MSR, as the metastable condensate 
will tunnel through the intermediate barrier. The macro- 
scopic tunneling rate is calculated semiclassically by the 
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WKB tunneling formula: 



-^tunnel 
N 



iV!/exp(-2 / y / 2[oj Q (r) - E] dr) 



= Nvexp(-2a) 



(16) 



where the limits of integration r\ and r2 are the inner 
and outer turning points of the intermediate barrier on 
the right (RIB) of oJo(r), E is the energy of the metastable 
condensate and exp(— 2a) is the WKB tunneling proba- 
bility. The frequency of impact (y) of the condensate on 
the RIB is approximately given by 



dr 



y/2 [E-Lu (r)} 



(17) 



where tq and r\ are the classical turning points of the 
metastable region. As TV increases, the net effect of the 
negative anharmonicity increases fairly rapidly. Hence, 
even though the minimum and stiffness of u>o(r) in- 
creases with TV, the difference (Aw) of the maximum 
of RIB (uJmax) and the minumum of the MSR (u> m i n ) 
decreases with increasing TV. Consequently, RIB disap- 
pears (Aoj = 0) when TV exceeds a critical value, N£ rst 
(to distinguish the critical numbers associated with the 
right side and the left side (see later) barriers, we name 
them as N^ rst and N^ cond respectively). This causes a 
new type of instability and eventual collapse. The tun- 
neling rate is appreciable only when E is close to uj max , 
and it increases rapidly as E approaches uj max . In Fig. 2, 
the tunneling rate (r^ n ™ ei ) of the lowest resonance state 
is plotted against the number of condensate atoms close 
to the critical point for various values of anharmonic dis- 
tortion. Near the criticality TV ~ N£t rat , the macroscopic 
tunneling is quite high and observation of this tunneling 
may be possible experimentally. The sharp peak near the 
criticality is attributed to the fact that the energy of the 
resonance state is close to the barrier height. Note that 
with increasing anharmonicity the right side barrier be- 
comes lower, which makes T t ^ nnel larger. For example, 
for A = —1.75 x 10 -5 o.u., the lowest resonance state 
near the critical point has tunneling probability 30%, 
whereas for A = —2 x 10 -5 o.u. corresponding tunnel- 
ing probability increases to 74%. Consequently, TV/ r lrst 
decreases with increasing |A|. In Fig. 3, we plot the reso- 
nance wave functions for two values of TV close to critical 
point {viz. TV = 567 and TV = 568, the critical num- 
ber being TV/,? rst = 570). The wave function within the 
metastable region is large and it has a small oscillatory 
part just outside the RIB. This clearly signifies that a 
part of the wave function leaks. For better clarity, the 
rapidly oscillatory part of the wave function is shown 
magnified in the inset of Fig. 3. Note that the amplitude 
of the leaked part increases as TV increases. At the criti- 
cal point (Acj = 0) the metastable region disappears and 
the whole wave function leaks, which corresponds to the 
collapse. The picture is qualitatively same as observed 
for attractive BEC in pure harmonic trap. However the 
phenomena near the present collapse is a bit different 
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FIG. 2: (color online) Plot of r% nnel (in atoms per second) 
vs N for the lowest resonance state near the criticality for 
various values of A (in o.u.). 



from the commonly observed collapse of attractive BEC 
in harmonic trap. In the latter case the metastable region 
is associated with a deep attractive well on the left side 
of the MSR, the metastable condensate tunnels through 
the left intermediate barrier (LIB) near the origin and 
scttels down in the deep well to form clusters. In a typ- 
ical attractive condensate we have checked that the size 
of the well is ~ 0.05^m. Hence, due to the high two-body 
and three-body collision rates within such a narrow well, 
atoms form cluster. The width of this wavefunction in 
the narrow well is of the order of 0.003/im, which is the 
order of the size of the atomic cluster. But in the present 
case the atoms which escape by tunneling outward will 
form a non-condensed Bose gas. 
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FIG. 3: (color online) Plot of the resonance wave function (in 
o.u.) vs r (in o.u.) for TV=567 and TV=568 for A = -2 x 
10 _J o.u. and a 3C = .00433 o.u.. The oscillatory part of 
the wavefunction immediately outside the barrier is shown 
magnified in the inset. 

As a further study to observe transition from a reso- 
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nance state to a bound state, we calculate decay rates 
for different values of effective interaction |iVa sc |. By de- 
creasing the effective repulsive interaction, we effectively 
enhance attraction between the atoms. Table f clearly 
shows that even for very slow decrease in | Na sc \ , r t 1 ^ nnel 
decreases rapidly, and very soon reaches a vanishingly 
small value, which manifests the transition from reso- 
nance to a bound state. 



Table f : Decay rates of lowest resonance states for 
different A in a repulsive BEC (a sc — 0.00433 o.u.). 



A = -f x 10~ 5 o.u. 


A = -2 x 10~ 5 o.u. 


\Na sc \ 
(o.u.) 


j^tunnel 
1 N 

(atoms/sec) 


\Na sc \ 
(o.u.) 


r^tunnel 
1 N 

(atoms/sec) 


4.30402 


2202.6475 


2.45944 


9096.67862 


4.29969 


53.7088 


2.45511 


1738.7188 


4.29536 


L3083 


2.45078 


108.5416 


4.29f03 


0.1107 


2.44645 


1.5149 


4.28670 


0.0f58 


2.44212 


0.1794 



Our result is qualitatively similar to earlier findings of 
Moiseyev et.al where the transition from resonance 
to bound state was discussed. However the earlier 
calculations 0, H[ used singular delta function potential 
in the mean-field equation and a negative offset potential 
was required to facilitate conversion of a quasi-bound 
state into a bound state, in a three dimensional BEC. 
In our present calculation we need no such offset. This 
deviation from GP results is attributed to the use of 
a realistic interatomic interaction having a hard core 
repulsive part at shorter separation. Moreover, our 
results provide realistic aspects which are relevant to 
experiments. 



B. Attractive BEC 

The situation becomes more interesting for the attrac- 
tive BEC in a shallow trap. We choose a condensate of 
85 Rb atoms with a sc = —1.832 x 10~ 4 o.u. For a clear 
understanding, we plot the effective potential in Fig. 4. 
The intermediate MSR is now bounded by two neigh- 
bouring barriers, one on the left (LIB) and one on the 
right (RIB) of unequal height. On the left side of LIB, 
a deep and narrow attractive well (NAW) appears. In 
the same vertical scale, we could not plot this deep well; 
hence it is not shown in Fig. 4. The RIB is the effect of 
negative anharmonicity which basically corresponds to a 
finite optical trap, whereas the LIB is the effect of the 
negative a sc which is commonly seen for attractive BEC 
in pure harmonic trap. The heights of the two barriers 
very strongly depend on two factors: first the anharmonic 
parameter and second the effective attractive interaction. 



Basically there is a competition between these two ef- 
fects which causes the shape of the effective potential to 
change in a complicated fashion with the increase in N. 
So throughout our study we fix A = — 9.37x 10 -6 o.u. and 
the effective attractive interaction is tuned by changing 
the number of bosons. The metastable condensate will 
have a finite probability of macroscopic quantum tunnel- 
ing through both the barriers. 

We start with few hundreds 85 Rb atoms in the trap, the 
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FIG. 4: (color online) Plot of the effective potential UJoir) 
in o.u. against hyperradius r in o.u. for a condensate of 



5 Rb atoms with N = 2660, A = -9.37 x 10" 
sc = -1.832 x 1CT 4 o.u.. 



and 



two neighbouring barriers are quite high and tunneling 
of the condensate through either of them is negligible. 
We have checked that in such a situation there is no sub- 
stantial leakage of the condensate through the associated 
barriers. With further increase in particle number, we ob- 
serve that the metastable region gradually becomes flat- 
ter, the corresponding condensate wave function expands 
slowly. With this wave function we calculate the average 
size of the condensate (r av ) [HJ as the root mean square 
distance of individual atoms from the center of mass of 
the condensate and is given by: 



1/2 



1/2 



(18) 



where R is the center of mass coordinate. In fig. 5a, we 
plot r av as a function of N, for TV increasing from a few 
hundred to a few thousand bosons. We find that r av 
increases slowly as expected (as RIB decreases and LIB 
does not change substantially, and as a consequence, the 
wave function spreads outwards). Finally at N = 2460, 
RIB vanishes and there is no MSR to hold the conden- 
sate, we call it a partial collapse. This is also reflected 
in the sharp fall in r av at N = 2460 (Fig. 5a). Thus 
N£j rst = 2460. The associated tunneling rate r$ nne! 
near the first collapse is shown in Fig. 6a. Near N^ rst , 
the condensate is associated with a large tunneling 
probability. Thus N^ rst is associated with first branch 
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of metastable condensate. 
With further increase in particle number we observe 
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FIG. 5: (color online) Plot of average size of the attractive 
condensate r av (in o.u.) as a function of TV in the the anhar- 
monic trap (A = -9.37 x 10" 6 o.u. and a sc = -1.832 x 1CT 4 
o.u.) near the first [panel (a)] and second [panel (b)] critical- 
ity. 

that the MSR reappears at TV = 2605, the second branch 
starts to develop and LIB decreases gradually. This is 
due to the fact that the net attractive interaction now 
dominates over the effect arising from the anharmonic 
distortion, as the former increases as jV ( J ^~ 1 ) while 
the latter increases as TV. Due to substantial increase 
in attraction, both the height of LIB and the local 
minimum of C0o(r) decrease rapidly, compared with the 
decrease of the height of RIB. Hence the MSR revives. 
As LIB decreases, the metastable condensate shrinks 
inwards, and we observe its behavior quite similar to 
what is seen in a pure harmonic trap: r av decreases 
sharply with increasing TV, as seen in Fig. 5b. Unlike 
the first metastable branch, in the second branch the 
fall of r av is fairly sharp (note the difference in the 
horizontal scales in the two panels) and quicker collapse 
occurs at TV = 2667. Wc name this as the second 
criticality (N™ con ). We also observe that near the 
second critical point, the condensate wave function is 
associated with an oscillatory part in the left side of 
LIB. At TV > N^ cond , the entire condensate collapses 
into the deep well, forming clusters. The associated 



tunneling rate r f ^ nnel for the second metastable branch 
has been calculated using Eq. (16) and Eq. (17) with 
the limits of the integrations suitably changed and is 
shown in Fig. 6b. The physical explanation for the 
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FIG. 6: (color online) Plot of T^ nnel (in atoms/sec) against N 
near the first [panel (a)] and the second [panel (b)] criticality. 



appearance of two distinct mctastabilities is as follows. 
The attractive interaction lowers the effective potential 
(u!o(r)) and the amount of lowering increases with TV as 

— — . This lowering decreases rapidly as r increases. 
On the other hand, the anharmonic term also lowers 
the effective potential, but the corresponding lowering 
is appreciable only for large r and it increases with 
r as well as TV. The increase with TV being roughly 
proportional to TV. Hence, when |A| is not too small and 
TV increases from a small value (< N^ rst ), the lowering 
due to anharmonicity (first lowering) at large r is much 
larger than that due to the interaction (second lowering) . 
Hence, with increasing TV, u {r) decreases for large r, 
giving rise to the first metastability and the appearance 
of the right intermediate barrier (RIB). As TV increases 
RIB decreases leading to the first criticality with the 
partial collapse at TV = TV/ r irsi . With further increase 
of TV, the second lowering at a smaller r increases 
faster than the first lowering. This causes reappearance 
of the MSR, whose minimum now gradually moves 
inwards. As TV increases even further, the second 
lowering for smaller r increases very rapidly, inducing a 
deep attractive well and an intermediate barrier on the 
left (LIB) between this well and the MSR. The second 
lowering decreases very rapidly with increasing r and 
is not strong enough at the position of RIB to alter it 
appreciably. As TV increases further, LIB disappears 
and the second criticality with collapse at TV = j\f^ cond 
results. However the two branches are discontinuous 
in the range 2461 < TV < 2604 for the present choice 
of parameter sets. As we have said earlier, there is a 
competition between the two controllable parameters, 
viz. interaction and anharmonicity. Thus the existancc 
of the discontinuous metastable branch will strongly 
depend on the choice of interaction and anharmonic 
distortion parameters. Our present study considers only 
a particular value of anharmonic distortion. So further 
study with other values of A is essential. 



IV. CONCLUSION 

In summary, we have applied a correlated many body 
method in three dimensions, incorporating a realistic 
interatomic interaction (van der Waals potential) to 
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study metastable condensates confined in a trapping 
potential with a finite barrier. The potential is taken 
as a sum of a quadratic plus a quartic term, which 
approximates an optical dipole trap i.e. a harmonic 
confinement combined with a Gaussian envelope. We 
obtain the complete quantitative description of the decay 
process of the quasi-stationary condensate with both 
repulsive and attractive interatomic interactions. Due 
to the use of a realistic interatomic interaction together 
with interatomic correlations in the many-body wave 
function and consideration of a finite number of atoms 
in the trap, our results exhibit more realistic picture. 
For a repulsive BEC, the quasi-stationary condensate 
can be stabilized by controlling the effective two-body 
interaction (through a sc ) and also the anharmonicity 
of the trap. By employing the WKB approximation, 
we calculate decay rates of such systems, which would 
be possible to measure experimentally. However, in 
contrast with earlier findings, in our present calculation 
no offset potential is required for the transition from 
a quasi-bound resonance state to a bound state. This 
difference is attributed to the use of a realistic inter- 
atomic interaction having a hard core at short range 
which prevents a catastrophic singularity at the origin 
as in the GP theory and produces a deep but finite well 
on the left of the left intermediate barrier. On the other 
hand, for an attractive BEC, in addition to the usual 
metastable condensate in a harmonic trap, we observe a 
new metastable branch which appears only for an inter- 
mediate range of particle number below the critical value 
for collapse due to attractive interaction only. The new 
metastable branch is also associated with an eventual 
collapse, for which the critical number is N^ rst . We also 
determine the decay rates of the metastable BEC due 



to quantum tunneling from both the metastable regions. 
However, the transition between these two branches is 
discontinuous. We have observed that this discontinuity 
strongly depends on the distortion parameter. However 
as the decay rate of 85 Rb atom in the new metastable 
region is quite high, the experimental study of this new 
phase may be difficult. But this technical difficulty may 
be circumvented by the proper choice of the parameters. 

Prediction of two branches of criticality , in particular, 
the fact that the criticality associated with the right 
side barrier appears and then disappears as iV increases 
from a small value upto N^ rst and then beyond, are 
the most significant new physics outcome of this work. 
From a technical point of view, the use of a many-body 
theory, incorporating all two-body correlations in the 
many-body wave function and a realistic finite range 
interatomic interaction with a strong short-range repul- 
sion are the realistic features. Deviations from earlier 
results are attributed to these. 
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